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Abstract 

The HZETRN code has been identified by NASA for engineering design in the next phase of space exploration 
highlighting a return to the Moon in preparation for a Mars mission. In response, a new series of algorithms 
beginning with 2005 HZETRN, will be issued by correcting some prior limitations and improving control of 
propagated errors along with established code verification processes. Code validation processes will use 
new/improved low Earth orbit (LEO) environmental models with a recently improved International Space Station 
(ISS) shield model to validate computational models and procedures using measured data aboard ISS. These 
validated models will provide a basis for flight-testing the designs of future space vehicles and systems of the 
Constellation program in the LEO environment. 


1. Introduction 

Improved spacecraft shield design requires early entry of radiation constraints into the design process to 
maximize performance and minimize costs. As a result, we have been investigating computational procedures to 

allow shield analysis starting with preliminary design concepts through high-fidelity final design models (Wilson et 

al. 2003). Of particular importance is the need to implement probabilistic models to account for design 

uncertainties (Wilson et al. 2004) in the context of optimal design processes (Qualls et al. 2003). These 
requirements need supporting tools with high computational efficiency to enable appropriate design methods. Only 
the HZETRN code has so far been identified for this purpose. As a result, Wilson et al. (2005) have prepared a 
review of past HZETRN code development, verification, and validation. Although there has been sporadic research 
to generalize this code over the last ten years, evaluation of code status among principal users demonstrated drift in 
the various code versions and databases. As a result of this renewed interest in HZETRN for future space systems 
design, there is now a systematic effort to advance, verify, and validate these codes in preparation to integrating 
them into engineering design processes. The present paper will describe the first several months of this renewed 
systematic effort. 

As NASA’s newly defined technology development spirals are now progressing, there is a need to identify 

a suitable code for the early spiral processes (development of a Crew Exploration Vehicle, CEV). We have chosen 
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the 1995 HZETRN version for this early spiral development. This code has had wide distribution and is largely the 
code used in the prior verification and validation processes described by Wilson et al. (2005). However, an 
assessment of the 1995 HZETRN code status demands the re-evaluation of the computational procedures and 
databases to account for software and database drift over the last ten years and the inherent limitations of these prior 
computational procedures. In the present paper, we review the basic derivations to identify oversights and improve 
solution stability and accuracy. In prior papers, we used both analytical and Monte Carlo benchmarks to both 
identify problems among various versions. We will distribute these benchmarks with the code package to evaluate 
future compiler and platform dependent problems in addition to future drift for whatever reason. This revised 1995 
HZETRN code, chronicled as the 2005 HZETRN (Wilson et al. 2006a), is recommended for standard practice and 
the present paper will focus on using measurements aboard ISS and recently improved models of the low Earth 
orbit (LEO) environment (Wilson et al. 2006b) for code validation. In addition to adding directional dependence 
and solar modulation to the improved trapped proton model AP8, we further developed a directional dependent 
geomagnetic transmission model. The main effect of this directional dependence of geomagnetic transmission is 
the penetration of particles below the vertical cutoff transmission model used in the past. Furthermore, we provide 
for the dynamic behavior of the Earth’s magnetic field covering the period 1945 to 2020. The emphasis of the 
present paper is to validate the modulated environmental models and 2005 HZETRN shielding code using detailed 
measurements on ISS. 

2. 2005 HZETRN 

The specification of the interior environment within a spacecraft and evaluation of the effects on the 
astronaut is at the heart of the space radiation protection problem. The Langley Research Center has been 
developing such techniques and an in-depth presentation is given in Wilson et al. (1991), although considerable 
progress has been made since that publication. The relevant transport equations are the coupled linear Boltzmann 
equations for a closed convex domain (Fig. 1) derived on the basis of conservation principles for the flux density 
(particles/cm 2 -sr-s-A-MeV) <p.(x, Q, E) for particle type j as 

Q»V<p(x,Q,E)= Zjo<Q,Q',E,E') <p,(x,Q',E') dQ'dE' - o(E) <p(x,O t E) (1) 

J Jk k J J 

where o(E) and o (Q,Q',E,E') are the shield media macroscopic cross sections. The a.(Q,Q',E,E') represent all 

j jk jk 

those processes by which type k particles moving in direction Q' with energy E' produce a type j particle in 
direction Q with energy E (including decay processes). Note that there may be several reactions that produce a 
particular product, and the appropriate cross sections for 
equation (1) are the inclusive ones. Exclusive processes are 
functions of the particle fields and may be included once the 
particle fields are known. Note, we will at times loosely refer 
to </>(. x, Q. E) as either flux or fluence and the usage should be 
clear from the context. The time scale of the processes in 
equation (1) are at most on the order of microseconds while 
time scales of boundary conditions are on the order of minutes 
or longer leaving the resulting interior fields in equilibrium 
with the particles at the bounding surface. 

The total cross section o(E) with the medium for each 

j 

particle type is 

or (E) = o°< (E) + o: el (E) + a r (E) (2) 

where the first term refers to collision with atomic electrons, 
the second term is for elastic scattering on the nucleus, and the 
third term describes nuclear reactions where we have ignored 
the minor nuclear inelastic processes (excited single particle 
states) except for low energy neutron collisions. The 
corresponding differential cross sections are similarly ordered. Many atomic collisions (~10 6 ) occur in a centimeter 
of ordinary matter, whereas ~10 3 nuclear coulomb elastic collisions occur per centimeter, while nuclear scattering 
and reactive collisions are separated by a fraction to many centimeters depending on energy and particle type. We 
include in the a '(E) term the nuclear decay processes. Solution methods first use physical perturbations based on 
the ordering of the cross sections with the frequent atomic interactions as the first physical perturbation with special 
methods used for neutrons for which atomic cross sections are zero. The first physical perturbation to be treated is 



Fig. 1 Geometric relations of quantities useful in solving 
equation (1). 
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the highly directed atomic collisions with mean free paths on the order of micrometers as observed in nuclear 
emulsion. The usual approximation is the continuous slowing down approximation leading to well specified range- 
energy relations but neglects the energy straggling that will be discussed later (see also Wilson et al. 2002a). The 
next term is the highly directed multiple coulomb scattering and is usually neglected in many models but is of great 
importance in understanding the transport of unidirectional ion beams leading to beam divergence and is treated in 
detail elsewhere (Shavers et al. 1993). The remaining nuclear reactive processes have been given main attention in 
past code developments. We will now consider the formal development of the relevant equations for further 
consideration. 


2.1 Continuous Slowing Down Approximation 


The collisions with atomic electrons preserve the identity of the ion and the differential cross sections are 

given as 

c r jk a>(Q,Q',E,E') = Z n cf‘. n (E’) 8(0 - O') 8 jk 8(E + Ar I e n - E’) (3) 

where n refers to the atomic/molecular excited states with excitation energies e n including the continuum. Note, the 
factor A- 1 results from the units of £ of A MeV (equivalent unit of MeV/nucleon). Although the atomic/molecular 
cross-sections cf.JE) are large (== 10 16 cm 2 ), the energy transfers e n are small (= 1-100 eV) compared to the particle 
energy. The atomic/molecular terms of equation (1) may be written as 


Zf<f l J 0,0', E,E') <p(x,Q',E') dO'dE'- cf‘.(E) <p.(x,Q,E) = £„ cf‘. (E + Aj'eJ <p(x,Q, E + A/'eJ - <f l .(E) <p(x,Q,E) 


■ Z n {(f‘ jn (E) <p.(x,Q, E) + A/e,, d E 

■ d E [S j (E)<P j (x,O,E)] + 0(e n 2 ) 


[cf‘ jn (E) <p(x,Q, E) ]} - (f. (E) <p.(x,Q,E) + 0(s„ 2 ) 


(4) 


where the stopping power S/E) is given as the sum of energy transfers and atomic excitation cross sections as 

S ] (E) = 2 n e n d u . n (E) (5) 

The higher order terms of equation (4) give rise to straggling and are neglected in the continuous slowing down 
approximation (csda) but are discussed elsewhere (Wilson et al. 2002a, Tweed et al. 2005). Evaluation of the 
stopping power by equation (5) is deceptively simple in that all of the excited states including continuum states of 
the atomic/molecular system need to be known. Furthermore, the projectile remains a bare ion except at low 
energies where the projectile ion atomic orbital states begin to resonate with the electrons of the media leading to 
electron capture and lowering of the ion charge. These issues are further discussed in Wilson et al. (1991) and Tai 
et al. (1997) in regard to the atomic data base. Equation (1) can be written in the csda as 

Q*V<p.(x,Q,E) - Aj l d E [Sj(E) 4>(xO, E) ] =Zfo Jk (Q,Q',E,E') <t> k (x,Q',E') dQ' dE' - a/E) <p.(x,Q,E) (6) 


where the right-hand side of equation (6) excludes the atomic/molecular processes now appearing on the left as an 
energy shifting operator in addition to the usual spatial drift term. Neutral particles would have null atomic cross 
sections for which the stopping term of equation (6) does not appear. Application of csda in both laboratory and 


space shielding has been wide and the resulting errors are 
discussed elsewhere (Tweed et al. 2005). Equation (6) can 
be rewritten as an integral equation (Wilson 1977) as 

<t>/x,a,E) = {Sj(E)Pj(E) <p.(r(Q,x),Q,E r ) 

+E/^dEA,P/E)fCf 4n dE"dQ' o jk (Q,Q',E',E') 
X<t> k (x+[Rj(E) - R j (E)]Q,Q’,E")}/S j (E)P J (E) (7) 

where T is the point on the boundary connected to x along 
direction -Q, E = Rj'[p - d + RJ, p is the projection of x onto 
Q (Fig. 1), d is the projection of r onto Q, Rj(E) is the 
distance an ion of type j of energy E will travel before losing 
all of its energy to excitation of atomic electrons, and Pj(E) 
is the probability a type j ion of energy E will have a nuclear 
reaction in coming to rest in the media. The usual range- 
energy relation is given by 

R/E) = fAj dE 7S/E ') (8) 



Fig. 2 Probability of nuclear reaction as a function of ion type 
and energy. 3 


(9) 


The nuclear attenuation function (Fig. 2) is given by 

P/E) = exp[- fAj cfj (E') dE'/S/E')] 

where the integral domains in equations (8) and (9) extend over the full energy range {0, E}. 

2.2 Straightahead Approximation 

The approach to a practical solution of equation (7) is to develop a progression of solutions from the simple 
to the complex allowing early implementation of high-performance computational procedures and establishing a 
converging sequence of approximations with established accuracy criteria and means of verification. The lowest 
order approximation using the straightahead approximation was guided by the nucleon transport studies of 
Alsmiller et al. (1968) using the Monte Carlo methods in which the differential cross sections were approximated as 

d jk (Q,Q',E,E ') - d jk (E,E ') 8(0 -O') (1 0) 

resulting in dose and dose equivalent per unit fluence to be within the statistical uncertainty of the Monte Carlo 
result obtained using the fully angle dependent cross sections. The relation of angular dependent cross sections to 
spacecraft geometry in space application was further examined by Wilson and Kandelwal (1974) using asymptotic 
expansions about angular divergence parameters demonstrating errors in the straightahead approximation to be on 
the order of the square of the ratio of distance of divergence to radius of curvature of the shield (a small error in 
most space systems). 

Equations (6) and (7) were examined for HZE ions using the following form for the projectile 
fragmentation cross sections as 

d Jk ( Q.Q'.E, E')- djk(E') N, exp{ - [QVe -OVe'] 2 /2 e jk 2 } ( 1 1 ) 

where d jk (E') is the cross section for producing fragment j from ion k, N, is the normalization constant for the 
exponential function, and s jk is the momentum dispersion parameter in the reaction. Substituting the interactive 
form of equation (11) into the integral term (target fragments treated separately, Wilson 1977) of the Boltzmann 
equation (6) yields 

ll d jk (Q,Q',E,E) <P k (x,Q',E) dQ' dE '= £ d Jk (E) {(j> k (x,Q,E) + Ed E <\> k (x,Q ,E) V[e jk 2 /(2mE)] 

+ Q»d a <p k (x, Q ,E) e jk /{ ( 2 mE) } (12) 

where the second term on the right hand side of equation (12) results from corrections in assuming the velocity of 
the ion is preserved in the interaction (Curtis and Wilkinson, 1972, Wilson 1977) and the third term is error 
resulting from the straightahead assumption (Alsmiller et al. 1968, Wilson 1977). The surprising result is that the 
velocity conserving assumption is inferior to the straightahead approximation for the nearly isotropic space 
radiation. There are great simplifications using equations (4) and (12) as given below 

Q*V(p.(x,Q,E) - A- 1 d E [S/E) <t>fx,Q, E)] = 1 d./E) <p/x,Q,E) - d/E) <p/x,Q,E) (13) 

which is strictly applicable to the HZE ions (Z > 2). The light ions and neutrons have additional complications 
arising from the broad energy spectra associated with their production, although the more favorable straightahead 
approximation is useful as indicated in equation (12) and as found by Alsmiller et al. (1967, 1968) in their 
verification process. The corresponding light ion (and neutron) Boltzmann equation is 

Q*V<t>/x,Q,E) - A/ d E [S/E) <p/x,Q, E)[ = If d./E, E') <p/x,Q,E') dE' - d/E) <p/x,Q,E) (14) 

where we have made use of the straightahead approximation as given by equation (10). Equations (13) and (14) 
have sufficient simplicity to allow an approach for both space and laboratory applications. The main force of the 
laboratory applications allow detailed model testing of the many atomic/molecular and nuclear processes. 

2.3 Marching Procedures 

The first version of the HZETRN code is based on the solution of equations (13) and (14) as guided by the 
Monte Carlo studies of the straightahead approximation by Alsmiller et al. (1968). We again specialize to solution 
along a ray Q directed along the x-axis for which equation (14) becomes 

[ d x - Af d E S/E) + d .(E)] f/x, E) = Ijd./E,E ') f(x,E ') dE' (15) 
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where d (E,E') are approximated for nucleons by the multiplicities of Ranft (1980), Bertini et al. (1972), and 
quasi-elastic contribution as described by Wilson et al. (1991). An immediate problem is the near singular nature 
of the differential operator, and transformation from energy to residual range coordinates as we did in developing 
the Green’s function greatly relieves this problem (Wilson et al. 1991). Unlike the Green’s function development, 
numerical procedures are simplified by introducing only a single residual range coordinate for all ions, and the 
residual proton range r is used as the common coordinate as 

r =f 0 E dE'IS(E') (16) 

and residual range of other particle types are related as v ; - r, ~ r with Vj = ZflAj which fails at low energies 
corresponding to low residual range due to electron capture into atomic orbitals characteristic to each ion type. The 
corresponding transport equation is 

[d x - Vj d r + 0 j(r)] ipj(x,r) = Zj r " (Vj /v k )s jk ( r ,r' ) % (x,r') dr' (17) 

where scaled flux is now (vjfor neutral particles such as neutrons are taken as unity in scaling relations, Wilson et 
al. 1991) 

yjj(x,r) = vj S( E) (j>fx, E) (18) 

and the scaled differential cross sections are 

S ik(r,r' ) = S(E) cf JE,E') (19) 

Errors in scaling of proton stopping and range parameters in arriving at the approximate transport equation (17) are 
compensated in part by solutions of equation (17) approaching a low energy equilibrium spectrum for ions given by 

Vj S(E) tj>.(x, E) => constant (20) 

where the constant is fixed by the higher ion energy. In distinction, the solution to equation (15) for ions has the 
low energy equilibrium spectrum 

A- 1 S(E) <p(x,E) => constant (21) 

also fixed by the higher energy flux for which the range scaling relation v, r, ~ r has better validity and the two 
constants are nearly equal so that equation (21) has improved accuracy over equation (20) at lower energies. This 
fact will require us to alter the flux unsealing relations as demanded by equation (21) to maintain accuracy at the 
lower energies. From equations (20) and (21), we can understand the simplicity of numerically solving equation 
(17) over numerical solution based on equation (15). The solution to equation (17) approaches a constant at small 
residual ranges allowing large separations in r grid values with smooth extrapolation to zero range whereas 
solutions of equation (15) vary as the nearly singular l/Sj(E) for which small E grid spacing is required leading to 
slow computational procedures. We tested the assumptions in equation (17) and unsealing according to relation 
(21) later in a separate report (Wilson et al. 2006a). 

The confusion caused by different scaling methods and associated coordinates for numerical procedures is 
justified by the simplification of the numerical representation of fluence of all particle types over a common 
residual range grid and simplification of the numerical procedures leading to high performance codes. Still a 
straightforward finite differencing of equation (17) can introduce unstable roots as had plagued the thermal 
transport problem for many years (Wilson et al. 1991). We will use the unconditionally stable methods of Wilson 
et al. (1991) arrived at by inverting the differential operator (Wilson et al. 1977, 1989, 1991) of equation (7) as 

ip/x.r) = exp[-£j(r,x)] ijj/0,r + Vj x) + ZJJ / + v/ v "exp[- £/ r,x ')]( v/v k ) s jk ( r +Vjx\ r ') ip k (x - x',r') dr'dx' (22) 

where the exponential is the integrating factor related to attenuation of the j type ions with 

U r,x) = f x Oj{ r+Vj X ') dx ' (23) 

and is related to equation (9). Equation (22) is a Volterra equation and can be solved either as a Neumann series 
(Wilson and Lamkin 1975) or with marching procedures (Wilson et al 1991). Note that the inverse mapping is 
taken as 

<t>j(x, E) = A. ipj(x,r)/S.(E) (24) 

to guarantee the equilibrium solution given as equation (21) at low energies away from the boundaries (note, the 
proton stopping power is used in case of unsealing the neutron flux). We assume in this treatment that the 



equilibrium constant resulting from equation (22) and given in equation (20) differs little from condition (21) for 
which the inverse mapping of equation (24) is most accurate. These approximations were tested by Wilson et al. 
(2006a). 

Two tracks are taken in implementing a marching procedure for equation (22) depending on particle type as 
demanded by the character of the nuclear processes. The problem naturally divides into “light ions” which will 
refer to all ions with atomic mass of four or less including neutrons and high charge-energy (HZE) ions having 
atomic mass greater than 4. The distinction arises from the energy and angle distributions of the double differential 
cross sections for which the HZE ions leaving a projectile fragmentation event have velocity nearly equal to that of 
the projectile as approximated by equation (1 1). Although the light ions are assumed to travel in the same direction 
as the projectile (see equation 10), they cover a broad energy distribution that cannot be ignored. The marching 
procedure is obtained by first considering equation (22) evaluated at x + h where h is the step size. Clearly, 

ip/x+hj) = exp[-£j(r,h)] %p j (x,r+v j h) + Ij 0 h f r + vjx “exp[-£ j (r,x')](v/v k )s jk ( r+ Vjx'y ) ipjx+h -x',r') dr'dx' (25) 

which may be used to develop a marching step from x to x + h once a means is found to approximate the field 
function ip/x,r) across the subinterval (x, x+h}. If h is sufficiently small such that a/r) h« 1 then, following 
lowest order perturbation theory (Wilson and Lamkin 1975), one has 

ip k (x +h - x',r') - exp[-E) k (r',h - x')] ip k [x, r' + v k ( h-x')] +0(h) (26) 

which may be used to approximate the integral in equation (25) giving results for the fields 0(h 2 ) as required to 
control the propagated error (Wilson et al. 1991). Substituting equation (26) into (25) and evaluating the 
attenuation factors at the interval midpoint (mean value theorem) results in 

ip/x+h,r) = exp[-Cj(r,h)] ip/x,r+Vj h) + Z k exp[-£/rJi/2) -£ k (r',h/2)] f r + yih/2 °° (v/vj F A jk (h,r,r) ip k (x,r'+v k h/2)]dr' 
+ 0(h 2 ) (27) 

where the integrand has been simplified using 

F A jk (h,r,r ) = f 0 h s jk ( r +Vjx' >r ') dx ' = F jk (r + Vj h, r ) - F jk (r,r) (28) 

and F jk (r,r) = f 0 c(r) dJE",E') dE" (29) 

with e(r) the energy associated with proton residual range r, and E' = e(r'). Note that if j corresponds to a neutral 
particle such as the neutron (/' = n ) then the above expressions are evaluated in the limit as Vj approaches zero in the 
range scaling relations resulting in the following (whereas the flux scaling factor for neutrons assumes v„ = 7) 

ip n (x + h,r) = exp[-o„(r) h] ip n (x,r) + Z k „exp[-o„(r) h/2 -£ k (r',h/2)] hf"(l/v k )sj r y ) r '+Vk h/2) dr' 

+ exp[- o„( r) h/2 - o n ( r ')h/2] h f " sj r y ) %p„(x, r ') dr ' (30) 

and similarly for the neutral k term (k - n ) when the j particle is charged 

ipj(x+h,r) = exp[-£j(r,h)]ipj(x,r+Vjh) + Z k „„ exp[ - r, h/2 ) - £ k ( r j h/2 )]/'( Vj/vi< ) F A jk (h, r, r '+ Vjh/2 ) 

xip k [x , r'+(vj+v k )h/2)] dr' + exp [-^/r, h/2) -o„(r') h/2] f r “’v i F A jn (h,r,r'+Vjh/2) ip n (x, r'+Vjh/2) dr'(3 1) 

where v„ in the flux scaling relation (24) is taken as unity. Equations (30) and (31) are solved on an equal-spaced x- 
grid Ax = h apart and a logarithmic spaced r-grid on two subintervals. The remaining integrals in these equations 
are approximated by (Wilson et al 1991) 

Jr k °° K(r n ,r) ipj(x,r') dr' » Z l=k K[r n ,(ri+ri +I )/2] f r [ ,+I ip/xp') dr' (32) 

where oo denotes a chosen upper limit tailored to the specific boundary condition. Note that the matrix of k"- values 
can be evaluated once on the r-grid and stored for subsequent steps providing high computational efficiency. 
Equations (30) and (31) provide the basis for the light ion transport of the HZETRN code. The HZE ion projectile 
( Aj > 4) coupling to the light fragments is contained in equations (27) to (31). 

The HZE fragments are produced with nearly the same velocity as the projectile ion as expressed in 
equation (13) and results in the simplified Boltzmann equation as 

[ d x - A;‘d E S/E) + d/E)] <p/x, E) = y k d./E) <p/x,E) (33) 

for which the scaled equations result in contributions from all HZE ions (with A k > 4) as 
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(34) 


ip/x,r) = exp[-£j(r,x) ] ipj(0,r + Vjx) + Ij 0 x exp[-^(r,x')] (Vj/v k ) cfjr+Vj x) ip k (x -x',r +Vjx') dx ' 

The corresponding marching equation (Wilson et al. 1991, Shinn et al. 1992) is given as 

ipj(x+h,r) = exp[-£j(r,h)] ipj(x,r+Vjh) + Ij () h exp[-^( r,x ')]( Vj/v k ) o jk ( r +Vjx) ip k (x+h-x', r +Vjx') dx ' (35) 

for which the integrand can be approximated for sufficiently small h using 

ip k (x +h - x',r+Vjx') = exp[-^ k (r+Vjx',h - x')] ip k [x, r+Vjx' + v k (h-x') ] +0(h) (36) 

allowing the following simplification 

ipj(x+h,r) = exp[-Cj(r,h)] ipj(x,r+Vjh) + Ij 0 h exp[-^(r,x')-^ k (r +VjX',h - x')](v/v k ) 

X o jk ( r+ Vjx) ip k [x, r+ VjX ' + v k ( h-x ')] dx ' (37) 

To evaluate equation (37), we use the mean value theorem that guarantees linear terms of the final integral to be 
zero. First, we expand the attenuation factor as 

L,j(r,x') = ff a/r+V/x") dx" ~ f x [o/r+Vj h/2) + d r Oj(r + Vj h/2 ) Vj (x"- h/2)] dx" (38) 

and similarly for 

tk( r+Vjx \h - x ') = J h - X ' a k [r+ Vj x '+ v k ( h-x ")] dx " ~f 0 h x ' [ a/ r+ Vj x '+ v k h/2) 

+ 8,.of r+Vjx '+ v k h/2 ) v k (x "- h/2 ) ] dx " (39) 

while applying the mean value theorem to the remaining factors of equation (37) and neglecting all but linear 
expansion terms in the integrand yields 

ipj(x+h,r) = exp[-tj(r,h)] ip/xj'+Vjh) + I k (v/v k ) o jk ( r+ Vjh/2) ip k [x, r + (Vj + v k )h/2] 
xfo exp{-Oj( r+Vjh/2 )x '-o k [ r+(Vj+ v k )h/2 )( h-x ')]} dx ' 

= exp[-C/rJi)] ipj(x,r+Vjh) + I k (v/v k ) a jk ( r +Vjh/2) ip k [x, r+(V/ + v k )h/2] 

x [exp{- Oj( r+ Vj h/2)h}-exp{- a k [r+ ( Vj+ v k )h/2)]h}]/{ o k [r+( Vj+ v k )h/2)] - Oj( r + Vj h/2)} + 0(h 2 ) (40) 

to be compared with the original HZETRN algorithm to 0[(Vj-v k )h] derived by Shinn et al. (1992) given as 

ipj(x+h,r) ~exp[-£j(r,h)] ip j (x,r+v j h)+I k (v/v k )o jk ( r ) ip k (x,r+Vjh) {exp[-Oj(r)h]-exp[-a k (r)h]} i[o k (r) - Oj(r)] (41) 

In earlier versions of BRYNTRN for proton/neutron transport, the flux scaling relation was taken correctly as 

xp j (x,r)=S(E)<P{x,E) (42) 

but carried over to the final BRYNTRN for light-ions/neutron transport in the most recent version of BRYNTRN 
(Cucinotta 1993). In coupling to HZETRN with scaling given by 

ipj(x,r) = vj S( E) <p.(x, E) (43) 

there is an inconsistency in flux scaling which must be accounted for. The appropriate coupling is given in 
equations (37) through (41) with the added factor of Vj/vk in the field coupling terms. The main effects on solution 
of the Boltzmann equation are expected for the light ions of H 2 , H 3 , and He 3 with only minor effects on the major 
light-ion/neutron components (n, H 1 , He 4 ). Note the main difference in the original and corrected code is 
determined by the ratios of Vj/vk of the field coupling terms. Current code verification processes are discussed in 
Wilson et al. (2006a). 

3. LEO Environmental Model 

The LEO environment consists of three main sources. Galactic cosmic rays (GCR) that penetrate the 
geomagnetic field, albedo neutrons from GCR interaction with the Earth's atmosphere, and particles trapped in the 
geomagnetic field. Three primary limitations in the environmental models are that the AP8 MIN&MAX are 
time/direction independent and that the vertical geomagnetic cutoff is used to describe the transmitted galactic 
cosmic rays. These models have been relatively successful in describing the radiation environment aboard the 
highly maneuverable Shuttle wherein anisotropies tend to be averaged. Such models will not be adequate for the 
ISS mainly oriented in the local horizontal plane along the velocity vector (minimum drag) except during battery 
charging. In a previous paper, we developed a dynamic/anisotropic trapped proton environment and general 
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geomagnetic cutoff model (Wilson et al. 2006b). These models are placed in a suitable form for evaluation of the 
incident radiation on the bounding surface of the six-degree of freedom motion of an orbiting spacecraft for shield 
evaluation. 

3.1 Albedo Neutrons 

Albedo neutrons result from the interaction of cosmic rays with the Earth’s atmosphere. As the cosmic ray 
intensities are modulated by solar activity so are the atmospheric neutrons modulated with time. The atmospheric 
neutron model is a parametric fit to data gathered by Langley Research Center studies of the radiations at SST 
altitudes in the years 1965 to 1971 covering the rise and decline of solar cycle 20. Scaling of the data with respect 
to geomagnetic cutoff, altitude, and modulation of the Deep River Neutron Monitor was found to allow mapping of 
the environment to all locations at all times resulting in an empirically based model for atmospheric neutrons 
(Wilson et al. 1991). The leakage flux F L (E) is closely related to the differential flux tf>(E,Q ) at the top of the 
atmosphere as follows 

F l (E) = f <t>(E,Q) cosddQ (44) 

where cos 6 is the direction cosine of the velocity vector with the zenith (note, (j)(E,Q) - 0 for cos 6 < 0). There are 
unresolved differences among various measurements of the leakage flux that is in part the assumed angular 
dependence of the differential flux. 


Table 1. Ratio of leakage flux to integrated flux as a function of energy and geomagnetic cutoff. 


Pc, GV 

Leakage to integrated flux ratio for energy (MeV) of- 


<.001 

.001 - 1 

1 - 10 

10- 19 

19 - 100 

100 - 2000 

0 

0.07 

0.63 

0.55 

0.48 

0.45 

0.38-0.23 

17 

0.07 

0.64 

0.56 

0.49 

0.53 

0.46-0.39 


The ratio of leakage flux to integrated flux is given by the New York University (NYU) group (Korff et al. 1979) in 
Table 1 that we use in the present model. The leakage flux was measured on the OGO-6 satellite by the University 
of New Hampshire (Lockwood 1972) during June 1969. The leakage flux spectrum was measured (Jenkins et al. 
1971, Preszler et al. 1972) and calculated (Korff et al. 1979) by various groups over Palestine, Texas during 
September 1971 and approximated by Wilson et al. (1999) as 


\0. 065/E for E £ 10 MeV 

F l (E) = \ (45) 

[0.0026 exp(- 0.01 IE) for E >10 MeV 

The total leakage flux is shown in Fig. 3 at various 
solar maxima and minima and in particular with the 
measurements on the OGO-6 satellite during June 
1969 which was a local maximum of cycle 20 (i.e., 
local minimum of DRNM). The leakage flux at the 
top of the atmosphere is extrapolated to low Earth 
orbit altitude according to Gauss’ law (varies as r" 2 ). 
There are some inconsistencies among various groups 
on specific neutron field related quantities. The 
present model is biased to the Palestine spectrum of 
September 1971, the OGO-6 latitude dependence 
during the local DRNM minimum of June 1969, and 
solar cycle dependence of the high altitude 
atmospheric neutron measurements as is clear from 
the above discussion (Wilson et al. 2002b). 



Fig. 3 GCR induced neutron leakage flux for different solar cycle 
maxima and minima compared to measurements by Univ. New 
Hampshire (Lockwood et al. 1972) 


3.2 Trapped Proton Environment 

The trapped proton population is modeled as AP8 for 1965 solar minimum and 1970 solar maximum (Vette 
1991). These inner zone particles are injected through the decay of albedo neutrons as they leak from the Earth’s 
atmosphere into the magnetic bottle defined by the Stormer forbidden zones (Wilson 2000). The inner zone 
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particles are lost from the trapping region by interaction with the tenuous atmosphere and generally have long 
trapping lifetimes above the fringes of the atmosphere. The inner zone consists of both proton and electron decay 
products. The average kinetic energy of the inner zone electrons is a few hundred keV. The electrons are easily 
removed from the spacecraft interior by the slightest amount of shielding and are mainly of concern to an astronaut 
in a spacesuit and materials external to the spacecraft. Within any pressure vessel such as has the Shuttle or ISS, 
the electrons are easily shielded by the micrometeoroid/debris bumper and pressure vessel wall. Of the trapped 
particles, only the protons with energies near or above one hundred MeV are of concern to the interior environment 
of the Shuttle or ISS. 

The particles trapped in the geomagnetic field were modeled from data obtained during two epochs of solar 
cycle 20 (solar minimum of 1965 and “solar maximum” of 1970) and are used with the geomagnetic fields on 
which the B/L maps were prepared (McCormack 1988). The 1965 analysis using the magnetic field model of 
Jensen and Cain (1962) resulted in particle population maps AP8 MIN (Sawyer and Vette 1976). The 1970 
analysis using the magnetic field model GSFC 12/66 (Cain et al. 1967) extrapolated to 1970 resulted in the particle 
population maps of AP8 MAX (Sawyer and Vette 1976). These models are considered the best available global 
representations of the trapped proton environment. This includes the known uncertainties in the AP8 models of a 
factor of 2 in LEO applications. 

It was believed at one time that better estimates of particle environments could be gained by evaluating the 
population maps defined on invariant Mcllwain coordinates over current magnetic field conditions. This 
interpolation would, for example, contain the westward drift of the South Atlantic Anomaly (SAA) observed in 
recent years by Badhwar et al. (1996). However, it was recognized by the Shuttle dosimetry group (Atwell et al. 
1989) that large errors resulted from such a procedure that does not account accurately the altitude shifts with 
changing magnetic fields. It was concluded that the use of the particle population maps interpolated over the 
magnetic field model for which the population map was derived would provide the best estimates of the long-term 
orbital averaged particle environments even though the westward drift is not represented. The westward drift is 
often introduced as a rotation of geographic coordinates (0.3 degree/yr) without regard to modifying the magnetic 
field over which the AP8 models were derived (Heynderickx 1996). 

In a prior report (Wilson et al. 2006b), we used the conventional westward drift of 0.3 degrees with an 
initial realignment of location such that AP8 MIN and MAX would be centered in 1965 and drift together 0.3 
degrees per year into the future with an improved modulation function. A renewed analysis of drifting of fields will 
be addressed in a companion paper wherein the spatial distribution of the radiation is addressed. The present paper 
will address long-term averages in which field drift is unimportant (Atwell et al. 1989) and focus on re-evaluation 
of the modulation factor and angular dependence. 

The proton omni-directional flux spectrum at any location and time f p (r,0,cp,E,t) is then extrapolated using 
the following functional form of Wilson et al. (1999, 2006b, see also Blanchard and Hess 1965) 


f p (r,6,(j),E,t) =f„ n Jr,0,<p,E) exp{-a [DRNMx F I0J - DRNMx F I07 I MIN ]} 


(46) 


where f p , mi „(r, 0, (j),E) is the proton flux at solar minimum shifted to 
time t and a p is evaluated using the solar maximum f ppmx (r,d,(j),E) 
related to AP8 MIN&MAX models but with the latitude and 
longitude tentatively shifted after centering in 1965 and 1970 as 
discussed by Wilson et al. (2007). In equation (46), the quantity 
(DRNMx F IO j) is averaged over the prior 14 months at solar 
minimum and 2 months at solar maximum with linear interpolation 
as determined by a best fit (Fig. 4) to the limited NOAAPRO model 
(Huston and Pfitzer 1998). Following the analysis of Wilson et al. 
(2007), we use the proton flux at solar minimum with 

fp,„u„(r,d,(j),E) = 0.5 f APSMIN (r,0 - 2.4 - 0.07At,cp + 4.1 + 0.19At,E) (47) 

and solar maximum with 

fp, ma .x( r ’@’ < P’E) = 0.6 f AP8MAX (r,0 - 4.8 - 0.07 At, cp +8+ 0.19At,E) (48) 





Fig. 4 Trapped proton environment in ISS 
orbit during June 2001. 


where At is time after the map epoch. The sunspot numbers and Deep River neutron monitor variation over solar 
cycle 23 (observed and projected) are used (Wilson et al. 2002b, Kim et al. 2006). 

Heckman and Nakano (1963) had studied the angular distribution of trapped protons with nuclear emulsion 
on rockets many years ago and presented a simple model of the pitch angle distributions about the geomagnetic 
field lines as related to the lifetimes of particles with guiding centers on different field lines. The protons’ velocity 
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vectors lie within 15 degrees of a plane perpendicular to the geomagnetic field line as seen in Fig. 5. Thus, those 
protons arriving from the east or the west differ in intensity according to the atmosphere scale height as related to 
the difference in population lifetime. The directional intensity is then defined for direction of arrival given by the 
local pitch angle 0 and azimuth X distribution as follows 


j(O.A) 


Ccxp 


-(jt/2-8) 


2o„ 


exp 


r cos /cos A 


(49) 


where r g is the particle gyro-radius about the field line, I is the dip angle, h s is the atmospheric scale height, a 0 is the 
width of the pitch angle distribution, and C is a normalization constant as discussed by Kem (1994). It was shown 
by Heckman and Nakano (1963) that a e depends on atmospheric scale height, altitude, and dip angle so that pitch 
angle distributions are nearly independent of particle energy. In distinction, the east-west asymmetry depends on 
the particle gyro-radius displaying marked energy dependence. The International Geophysical Reference Field 
(IGRF) is implemented at the current time and the scale height is found from the solar modulated fit of Pfitzer 
(1990) used by Badhwar (1999) to organize the Shuttle dosimetry data and given as 

p(r) = p 0 exp{-(h-120)/[A (h-103) 1/2 ]} (50) 

where po = 2.7 x 10' 11 g/cm 3 , h is altitude in km, and fitting parameter A = 0.99 + 0.518 [(F + Fbar)/1 10] 1/2 with 
Fbar the average of F given as the FI 0.7 solar radio output parameter over three prior solar rotations (81 days). The 
trapped protons are encountered by LEO spacecraft in the SAA, and ISS encounters this region from two directions 
(first with ascending node and second with descending node) as 
occurs during orbit precession. In addition, the radiation 
incident on the outer surface of the spacecraft is required for 
shield evaluation, and the attitude of the spacecraft is never 
fixed but limit cycles if not under reorientation from required 
maneuvers. The angular distribution averaged over spacecraft 
attitude in the region of radiation encounter is then to be 
evaluated. This is accomplished by relating the orientations in 
the spacecraft frame through yaw plus heading, pitch, and roll to 
the local vertical reference frame where the radiation 
environment is evaluated. We normally use 970 ray directions 
to evaluate the boundary conditions for shield evaluation and 
will use these same directions for evaluation of the directional 
environment. 

3.3 Geomagnetic Transmission Factor 

In the past, the geomagnetic transmission factor used 
was based on extrapolation of a world map of vertical cutoff rigidities of Smart and Shea (1983). In this model, it 
was assumed that there is no transmission below the vertical cutoff and 100 percent transmission (excepting the 
Earth’s shadow) above the vertical cutoff. There is, in fact, variable transmission dependent on angle of incidence 
relative to the east direction. The directional dependent cutoff for a centered dipole field is given by Stormer theory 
in terms of the cutoff rigidity as 

P C (Q) = M cos 4 (X)/{r 2 [1+(1 - cos 3 (X) cos(w E )] 2 } (51) 

where Q is the local direction of incidence, M is the dipole 
moment, X the magnetic latitude, cos((jo e ) is the projection of Q 
on the eastward direction. In application to the geomagnetic 
field, one must first transform the geocentric location into 
eccentric dipole coordinates. It was shown by Quenby and 
Webber (1959) that higher multipole contributions can be 
approximated by replacing X of equation (51) with X’ given by 

X’ = tan -1 {[V c + 0.52 5 V ]/[2(H C + 0.52 b H )]} (52) 

where V c and H c . are the vertical and horizontal field components 
at the location and 5 V and 8 H field deviations from dipole values. 

Fig. 6 Angular distribution of geomagnetjg 
cutoffs at the equator near SAA. 




Fig. 5 Angular distribution of 82 MeV protons 
near center of South Atlantic Anomaly. 


The Quenby/Webber model is of itself rather inaccurate (up to 25 percent) but a renormalization of the transmission 
factor (along the vertical cos(o) E ) = 0) using the more accurate vertical cutoffs as evaluated by Smart and Shea 
allows reasonably accurate anisotropic cutoff values (De Angelis et al. 2003). One result of the new transmission 
factors will be the admission of lower energy cosmic rays than in the currently used vertical cutoff model. In the 
present model, we use the IGRF field model evaluated for arbitrary dates from 1945 to 2020. A typical 
transmission factor of the model is shown in Fig. 6 at a location on the equator over the Earth near the SAA at 400 
km in March 1995. The western approaches are seen as the high cutoff regions on the sphere. The cosmic ray 
model is that of Badhwar-O’Neill (1995) extrapolated according to neutron monitor count rates (Wilson 2006b). 


4. Model Validation 


400 MeV/u Carbon beam 



Fig. 7 Liulin detector response model for incident 
400 A MeV carbon beam. 


Model validation has generally used Shuttle 
measurements with TLD100 to examine the long-term modulation 
of the trapped environment with more detailed studies with ISS in 
which the angular variations have some importance (Hugger et al. 

2003, Wilson et al. 2006b). The present step of model validation 
will focus on the ISS DOSMAP Experiment (E094, PI: Dr. 

Gunther Reitz, DLR) consisting of passive thermo-luminescent 
detectors (TLD), passive nuclear track detector packages 
(NTDPs), active dosimeter telescopes (DOSTELs), and mobile 
dosimeter units (MDUs, identified as Liulin-E094, PI: Dr. T. 

Dachev 2006). The Liulin energy loss spectra have been 
correlated to the NTDP derived liner enetgy transfer (LET) spectra 
by dividing the MDU energy loss spectra by the Liulin detector 
thickness of 0.306 mm Si and conversion to equivalent water. 

Unlike the TLD and NTDP data (Benton 2004), the advantage of 
the active detectors is their temporal variations used to separate 
SAA and GCR contributions. The Liulin energy loss spectral 

response functions can be accurately modeled as a function of incident angle, particle type, mean energy loss, and 
energy loss fluctuations due to straggling (Wilson et al. 2002). The response for C-12 ions at 400 A MeV is given 
in Lig. 7 and compares favorably with the measurements using the HIMAC in Chiba (Uchihori et al. 2002). 

The data were recorded aboard 

ISS 6 A during the late June - early July Tahle 2 • Descriptive parameters for LIULIN data sets. 

time frame of 2001 for several locations in 
the US Lab Module and Node 1 at 30- 
second time intervals. The positional 
trajectory data were supplied along with 
the time sequential measurements. Table 
2 denotes pertinent information for the 
data sets analyzed in the present work. 

The recently-developed NASA-Langley space radiation environment code (GEORAD), described 
elsewhere (Wilson et al. 2006b, 2007), has been utilized to simulate the environments for the orbiting ISS for the 


Designation 


Dates of 
Measurement 


Measurement 
Time Span, hr 


MDU 


Number of 
Data Records 


SET-1 

July 6-9,2001 

82.13 

2 

9856 

SET-2 

July 6-13,2001 

138.60 

1 

16633 

SET-3 

July 6-13,2001 

138.55 

3 

16625 

SET-4 

June 26-28,2001 

57.46 

1 

6893 



Fig. 8. ISS photo in 6A configuration and CAD-modeled Lab Module with detector locations. 
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time periods and the orbital tracks associated with each data set. During the measurement period, ISS was in its 
designated 6A development stage shown in Fig. 8 with the CAD model simulation. 

The CAD model utilizes the I-DEAS® commercial software package and consists of 512 individual 
components (including detector parts) specified by their spatial location and effective density (mass/volume). The 
established CAD model then enables the construction of a directional thickness distribution about each target point. 
Finally, with a given environment (trapped protons and OCR) and an associated shield distribution, the NASA- 
Langley transport code 2005 HZETRN is utilized to predict particle flux and dosimetric quantities at each target 
point indicated in the figure using the E094 site designations. 

4.1 Mission Segment Data 

Sample records of relevant data 
(formatted as ASCII values) are given in 
Table 3 as part of the information relating to 
sequential time points along the orbital track. 

Note that the exposures have been expressed 
in terms of a dose rate (it Gy/hr). The general 
nature of the high quality and consistency of 
the data is shown in Fig. 9 where the first two 
days of data set 2 are shown as measured dose 
rate as a function of time (corresponding to 
MDU-1 located at the TLD103 site in Fig. 8). 

The repetitive low-level ripple is 
indicative of the GCR exposures, while the 
sharp “spikes” represent passage through the 
SAA. In order to extract the SAA trapped 
proton dose from the underlying GCR, an 
analysis script was written that effectively filters the SAA “spikes” and defines the time interval for SAA passage. 
A sinusoidal function has been used to fit the smoothed low-level GCR for the time period and subtracted from the 
measured total to obtain the SAA exposure for a single pass. An illustration of this procedure is shown in Fig. 10, 
where the dose rate data corresponding to an SAA pass (first “spike” of Fig. 9) is plotted over the appropriate time 
interval. The fit for underlying GCR is also shown. 

MDU-1 Data (6-13 Jun 01) 


Table 3. Sample Data Records from L1ULIN Data Set 
Elapsed Time Dose Rate Latitude Longitude Altitude Heading 


(hr) pGy/hr (deg N) (deg E) (km) (deg) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

19.0500 

3.41 

-43.65 

294.01 

390.36 

65.14 

19.0583 

2.92 

-42.62 

296.15 

390.21 

63.49 

19.0667 

3.77 

-41.55 

298.22 

390.05 

61.66 

19.0750 

3.99 

-40.43 

300.21 

389.88 

59.93 

19.0833 

7.92 

-39.28 

302.14 

389.71 

58.24 

19.0917 

16.62 

-38.09 

303.99 

389.52 

56.56 

19.1000 

26.74 

-36.87 

305.79 

389.34 

55.02 

19.1083 

45.20 

-35.62 

307.52 

389.15 

53.35 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 
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Fig. 9. Sample data from LIULIN detector. 
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Fig. 10. Sample dose rate data for single pass through SAA with associated underlying GCR. 

The total mission segment dose is obtained by direct integration of the dose rate data with respect to time. 
The corresponding integrations of the individual SAA passes after GCR subtraction render the mission segment 
exposure for trapped protons. Finally, subtraction of the trapped contribution from the total segment dose results in 
the GCR dose. A summary of the measured data thus analyzed is given in Table 4. Note that each data set consists 
of SAA passes for which the ISS is either ascending (headed northward) or descending. Consistent differences for 
these two modes are observed, and result from the directional environment as well as location and directional 
dependent shield distributions. Such differences in dosimetry of the ascending and descending trajectories has been 
previously investigated by Dachev et al. (2006) and Hugger et al. (2003). 

4.2 Description of Model Calculations 

The positional data for each mission segment (time, latitude, longitude, altitude, heading plus yaw, pitch, 
roll) were used as input to the GEORAD environment code to calculate the fluence energy spectra for the GCR and 
trapped protons. For the calculations pertaining to this work, the cumulative directional fluence along 970 
directions has been evaluated in the vehicle coordinate system in the +XVV orientation: (x-axis along velocity 
vector, z-axis toward earth center, y-axis completing the right-hand system). Yaw, pitch and roll angles are 
assumed to be zero (neglects the small limit cycle motion), with heading direction specified as angle between 
direction of flight and geographic north (heading). Directional properties of the SAA protons have been accounted 
for throughout, and orbit-averaged directional transmission factors for GCR have been used. With respect to solar 
cycle conditions, the relevant data periods correspond very nearly to solar maximum conditions when GCR and 
trapped proton fluences are lowest. 

Results of fluence calculations integrated over all angles are presented in Figs. 11 and 12 for the trapped 
protons and GCR, respectively. The proton spectra have been plotted as a total fluence for each mission segment, 
and reflect the durations thereof. The GCR spectra are plotted in terms of a daily flux for several important GCR 
ions, and they illustrate the effect of flux reduction at lower energies due to geomagnetic cutoff. The flux spectra 
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are then used as input for the 2005 HZETRN transport code to generate dose vs. depth functions (see Figs. 13 and 
14). 


Table 4. Summary of Results from Data Set Measurements 


Data Set 1 

MDU#2 (6-9 July) 


Asc. 

Desc. 

4.5 

2.37 

15.98 

8.98 

0.66 

2.04 

1.35 

6.71 

14.65 

5.72 

5.73 

3.24 

7.98 

8.15 

14.05 

0.93 

2.85 


16.42 



SUM 84.17 38.14 
Final Sum 122.31 


Data Set 2 

MDU#1 (6-13 July) 


Asc. 

Desc. 

6.21 

4.86 

16.56 

11.52 

0.91 

2.24 

1.96 

2.38 

15.15 

9.27 

5.74 

6.59 

0.5 

1.35 

9.91 

10.25 

15.4 

0.96 

3.94 

1.07 

17.18 

6.78 

6.85 

7.97 

17.57 

5.34 

0.64 

10.71 

2.59 

1.99 

17.75 

2.26 

5.78 

1.9 


9.4 


6.02 


144.64 102.86 
247.5 


Data Set 3 

MDU#3 (6-13 July) 


Asc. 

Desc. 

4.64 

4.48 

15.98 

9.88 

0.96 

2.18 

1.28 

1.79 

14.2 

8.94 

6.31 

5.56 

8.31 

5.94 

14.41 

8.46 

2.67 

0.97 

16.4 

6.63 

4.76 

6.12 

15.96 

4.74 

0.64 

9.12 

1.42 

1.64 

14.47 

1.94 

5.57 

8.88 


4.83 


Data Set 4 

MDU#1 (26-28 June) 


Asc. 

Desc. 

2.75 

9.98 

15.22 

1.08 

1.97 

3.67 

0.84 

10.09 

11.43 

3.32 

8.49 

1.63 


8.06 


7.75 


127.98 92.1 
220.08 


40.7 


45.58 

86.28 


Total Mission Dose: 

444.13 693.2 709.14 284.41 

Total GCR Dose: 

321.82 445.7 489.06 198.13 


It is assumed that the ISS structure is primarily aluminum alloy (A1 2219), and particle flux transport is 
calculated for traverse through equivalent g/cm 2 A1 amounts. The resultant dose rates are computed for energy 
deposition in silicon (detector material) and shown in Figs. 13 and 14. It is seen that dose rates for the mission 
segments are very nearly the same for given aluminum thickness values for GCR with a small separation for 
trapped protons dependent on specific SAA crossings. For reference, the predicted free space solar maximum GCR 
dose rates are shown in Fig. 14 to display the effects of geomagnetic transmission factors. 

In order to calculate exposure at the specified target locations, directional distributions of shield amounts 
including the MDU materials are required. A spherical coordinate grid of 970 rays subtending equal solid angles 
has been implemented. The cumulative angle averaged shield thickness distributions based on the CAD-modeled 
ISS 6A configuration are shown in Fig. 15. Data measurements for the MDU#4 location were not available as of 
this writing, and the MDU#4 thickness distribution was not used in this study. When the dose vs. depth functions 
are used in conjunction with the directional thickness distributions, the total mission dose may be predicted. The 
results of these end-point calculations are given in Table 5; direct comparison with measurements is made in the 
following section. 
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Flux, nuclei/(cm A 2-AMeV-day) Fluence, protons/(cm A 2-MeV) 



Fig. 11. Calculated trapped proton fluence for data set simulation. 



Kinetic Energy, A-MeV 


Fig. 12. Calculated daily GCR flux for data period simulation. 
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Dose Rate, micro-Gy(Si)/hr 



Fig. 13. Transport calculation results (Dose Rate vs. Depth) for simulated data set trapped protons. 



Fig. 14. GCR transport calculations for simulated data sets. 
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Aluminum Amount, T, g/cm A 2 


Fig. 15. Cumulative angle averaged distributions of shield amounts for the detector locations. 


Another by-product of the 2005 HZETRN code is the 
calculation of linear energy transfer (LET) at the detector location. A 
commonly used representation of energy deposition per unit thickness 
for a given exposure condition is the integral LET spectrum that 
indicates particle flux resulting in a given energy deposited along the 
path through the detector. The combined (GCR + trapped protons) 

LET spectra for the four data sets are shown in Fig. 16. Although the 
DOSMAP data includes information sufficient to infer LET spectra, 
at the time of this writing the measured data were not readily 
available and remain a subject for future analysis. 

5. Comparison of Results and Discussion 

A summary chart of the direct comparison of exposure data measured with that calculated is given in Fig. 
17 in terms of the average hourly dose rate for each dataset. In general, agreement is fair to good, but despite all 
efforts to replicate the experiment, differences occur that suggest room for improvement. Since the measured data 
are of high quality and consistency, a first reaction is to examine possible weaknesses and uncertainties in the 
theoretical calculations. The trapped environment models used exhibit strong sensitivity to magnetic field 
specification and altitude (especially in LEO). The GCR component is a rather strong function of assumed solar 
cycle modulation between maximum and minimum values (in essence, also driven by geomagnetic field 
transmission factor). The comparison calculations are also strongly influenced by the modeled thickness 
distributions around the target points (detector locations), particularly with respect to the trapped proton exposures 
and further refinement of the ISS model may reduce differences. The distribution of shield thickness is also 
inherently coupled to the modeled directionality of the trapped protons in the SAA. The solar cycle modulation, 
shield amount distributions, and trapped proton directional flux have all been addressed in the present calculations 
with the best information available. 


Table 5. Final Total Mission Segment Dose 
Calculation Results, fxGy 

Data 

Set 

Trapped 

Protons 

GCR 

Total 

1 

193 

206 

399 

2 

240 

368 

608 

3 

104 

368 

472 

4 

111 

152 

263 
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Fig. 16. Particle flux LET spectra calculated for simulated data set environments and detector locations, 


I 


f 


□ Trapped 

□ GCR 



D-1 (m) D-1 (c) D-2 (m) D-2 (c) D-3 (m) D-3 (c) 

Data Set Identification 


D-4 (m) D-4 (c) 


Fig. 17. Comparison of measured (m) with calculated (c) average dose rates. 
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Another factor that could impact the calculations is that of the detectors’ directional response (field of 
view) which has been addressed in the present work as the asymmetrical distribution of the instrument mass but not 
with respect to the directional response of the silicon detector (see Fig. 7). This effect could be of significant 
importance for the less penetrating trapped proton component. These considerations suggest areas for examination 
in future analysis. It does appear that the GCR contributions of the four datasets as predicted by the models in Fig. 
17 are consistently low by about 15-20 percent and may be due to use of the older implementation Badhwar- 
O’Neill (1995) GCR model wherein the decades older data of the model are extrapolated using only current and 
historic Sunspot Numbers (Wilson et al. 2002, Kim et al. 2006), and the use of the new Badhwar-O’Neill GCR 
model (O’Neill 2006) is expected to improve the modeled result and reduce these differences. 

6. Concluding Remarks 

Consistent differences in the GCR comparison exist in that calculations result in lower predicted exposures 
for all mission segments. These exposures are sensitive to the strength of the solar maximum condition. A weaker 
solar maximum than that used in the present calculations would increase these predicted exposures. For the trapped 
exposures, differences are not consistent; very good agreement is indicated for data sets 2 and 4, whereas 
significant differences are seen for sets 1 and 3. For set 1, greater trapped exposure is predicted while for set 3 the 
trapped exposure is substantially underestimated. Since the trapped proton spectra are relatively “soft”, they exhibit 
high sensitivity to effective shield amounts. Additionally, the trapped environment is modeled with solar cycle 
modulation and imposed drift of the aging AP8MIN and AP8MAX static models. Future analysis of data such as 
that provided by the LIULIN system could well lead to up-grading the trapped environment with possible definition 
and inclusion of inherent short-term variability of flux intensities. 
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